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ABSTRACT 

When the accretion rate is more than a smaU fraction of Eddington, the inner regions of accretion 
disks around black holes are expected to be radiation-dominated. However, in the a-model, these 
regions are also expected to be thermally unstable. In this paper, we report two 3-d radiation MHD 
simulations of a vertically-stratified shearing box in which the ratio of radiation to gas pressure is ^ 10, 
and yet no thermal runaway occurs over a timespan ~ 40 cooling times. Where the time-averaged 
dissipation rate is greater than the critical dissipation rate that creates hydrostatic equilibrium by 
diffusive radiation flux, the time-averaged radiation flux is held to the critical value, with the ex- 
cess dissipated energy transported by radiative advection. Although the stress and total pressure 
are well-correlated as predicted by the a-model, we show that stress fluctuations precede pressure 
fluctuations, contrary to the usual supposition that the pressure controls the saturation level of the 
magnetic energy. This fact explains the thermal stability. Using a simple toy-model, we show that 
independently-generated magnetic fluctuations can drive radiation pressure fluctuations, creating a 
correlation between the two while maintaining thermal stability. 

Subject headings: accretion, accretion disks — instabilities — MHD — radiative transfer 



1. INTRODUCTION 



the midplane is 



As first demonstrated bv IShakura &: SunvaevI (|1973f ). 
radiation pressure must dominate gas pressure in the in- 
nermost parts of black hole accretion disks radiating at 
any non-negligible fraction of the Eddington limit. Al- 
though the original argument made use of the assump- 
tion that the stress and the pressure are related by a 
fixed ratio a, the critical radius within which radiation 
pressure exceeds gas pressure depends so weakly on the 
value of a that this conclusion would be virtually unal- 
tered even if the ratio of stress to pressure varied sub- 
stantially from place to place. In this sense, the ex- 
pectation of radiation dominance in the inner rings of 
brightly-radiating accretion disks rests only on the sup- 
position that stress and pressure are comparable, the di- 
mensional analysis foundation from which the a model 
began. It is therefore crucial that we understand the 
role that radiation pressure plays in the physics of ac- 
cretion disks if we are to understand luminous black hole 
sources. Yet, for the past thirty years the canonical inter- 
nal equilibrium of disks in the radi ation-dominated limit 
has been believed to be unstable (iLightman fc EardlevI 
[1971: IShakura fc SunvaevlfTOT^ lPiranlll97af) . 

In contrast to other physical regimes, the vertical struc- 
ture of a stationary disk is in some ways very tightly con- 
strained if radiation pressure truly dominates the hydro- 
static support of the disk against the tidal gravitational 
field of the hole. Neglecting relativistic corrections for 
simplicity, one finds that the outward radiation flux in a 
geometrically thin disk as a function of height z above 



Fiz) = -g{z) 



(1) 



where c is the speed of light, k is the opacity (usually 
dominated by electron scattering) , and the orbital angu- 
lar velocity Q determines the vertical tidal gravitational 
acceleration g. Because energy conservation in a time- 
steady disk demands (if we neglect inner boundary condi- 
tion effects) that the flux emerging from the disk surface 
is (3/87r)Mri^, the half-thickness H of the disk depends 
only on the accretion rate M and is independent of radius 
and the nature of the turbulent stress: H = 3kM / {8ttc) . 

Moreover, equation ([T]) combined with the condition of 
radiative equilibrium immediately implies that the dissi- 
pation rate per unit volume Q is constant as a function 
of height above the disk midplane: Q = dF/dz — j k. 
Given that Q must, on average, be given by the turbu- 
lent stress times the rate of strain r|c?ri/dr|, we also 
have a tight constraint on the vertically-averaged stress: 

2cO , . 

Tr, = (2) 

a conclusion first reached in slightly different form by 
IShakura fc SunvaevI (|1976[) . 

On the other hand, other aspects of radiation- 
dominated disks are left totally unconstrained. Because 
K (when it is dominated by electron scattering) is inde- 
pendent of density, the vertical density profile p{z) is 
completely irrelevant to both the hydrostatic and the 
thermal equilibrium of the disk. In the absence of ad- 
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ditional information about how the dissipation rate de- 
pends on density, it is therefore completely undeter- 
mined, and even within Shakura-Sunyaev models it must 
be specified in some ad hoc manner. Often it is assumed 
that the dissipation rate per unit mass is constant, which 
leads to a vertically constant density out to the photo- 
spheres. This element of arbitrariness is in addition, of 
course, to the ad hoc "a-viscosity" prescription for the 
stress used in such models. 

The freedom to choose arbitrary values of a in that 
stress prescription has always been one of its unfortu- 
nate aspects, but in the radiation pressure dominated 
regime we also have freedom in the choice of which 
pressure with which to sc ale the stress. The standard 
IShakura fc SunvaevI ()197 3^ assumption is that the stress 
should be proportional to the total (gas plus radia- 
tion) pressure, and this ch oice leads to both thermal (e.g. 
IShakura fc Sunvaevlll976l) and "visco us" (more properly, 
"inflow" ; iLightman fc EardlevI [l97l instabilities where 
the radiation to gas pressure ratio exceeds 3/2. The ther- 
mal instabilities generally have faster growth rates than 
the inflow modes, and are the ones of primary interest 
here. 

The origin of the thermal instability is most simply ex- 
pressed in terms of the different dependences of the heat- 
ing and cooling rates on the midplane te mperature T , as- 
suming constant surface mass density E (|Piranlll978f ). In 
the limit of radial wavelength long compared to the disk 
thickness, radiative diffusion implies that the cooling rate 
per unit area P- ~ AcaT^/{3KY,) cx T"*/!], where a is the 
radiation density constant. This must balance the heat- 
ing rate per unit area F"*" ^ 2HTrci,r\d^l/dr\. Combin- 
ing these two expressions with the a-prescription = 
and the condition of hydrostatic equilibrium, we 
find H = 2aT4/(3Sl}2) ^nd F+ ~ 2aa'^T^ / (Sni:) oc 
r^/S. Thermal instability follows because positive tem- 
perature perturbations lead to greater increases in heat- 
ing than in cooling. It is also possible to express these 
relationships in terms of the total radiation energy con- 
tent per unit area U ~ 2aT'^H, in which case F~ oc 
U^l"^ jY}!'^ and cx C/; once again, instability is indi- 
cated. 

However, it has never been clear whether the assump- 
tion that the stress is proportional to the radiation pres- 
sure is correct, and the existence of the thermal insta- 
bility has therefore always been suspect. Alternative 
prescriptions that make the stress proportional to ei- 
ther the gas pressure alone or some combination of the 
gas and radiation pressure have been advocated over 
the y ears, and give rise to more stable con figurations 
(e.g. ISakimo to fc Coronitil [l98ll: IStella fc Rosner 198 
[Rurm 1985; Sz uszkiewiczlll990t iMerloni fc Fabia'SlloO 
[Merloni, 2003, ). Moreover, all of these prescriptions re- 
fer to vertically-integrated or one-zone vertical structure 
models of the accretion disk. The actual vertical distribu- 
tion of dissipation may also be important. If much of the 
local accretion power is actually dissipated in optically 
thin regions above and below the disk, then the disk itself 
can become supported by gas pre ssure and there would 
then be no thermal instability (jSvensson fc Zdziarskil 
[1991 . 

It is now widely believed that the turbulence in black 
hole accretion disks resembles that seen in simulations 



of the nonlinear develo pment of the magnetoro tational 
instability (MRI, e.g. iBalbus fc Hawlevl Il998h . Such 
simulations, if they include radiation physics, could in 
principle help resolve these uncertainties by calculating 
the spatial and temporal structure of the turbulent dis- 
sipation within the disk, and checking to see whether a 
long-lived, and therefore presumably stable, equilibrium 
exists in the radiation pressure dominated regime. It 
may also be possible to investigate the scaling of aver- 
age turbulent stresses with local averages of gas and/or 
radiation pressure using the data from these simulations. 

Well-resolved, global simulations of optically thick 
accretion disks that include radiation transport have 
not yet been published, but it is computationally 
feasible to investigate the thermal physics locally 
within the disk u sing a stratifie d shearing b o x ge- 
ometry JF randcnb ure et all 119951 : IStone et all 119961 : 
[Miller fc Stone ,2000l l The necessity of (sheared)- 
periodic radial boundary conditions zeroes any net accre- 
tion rate through the box, precluding any investigation of 
inflow instabilities, but thermal instabilities can still be 
studied. These same boundary conditions also preclude 
studying any modes with wavelengths longer than twice 
the box width, except for infinite wavelength modes — 
the sheared periodic radial boundary conditions place no 
restriction on modes that are completely independent of 
radius. 

The first attempt at this was made by [Turneil 
(|2004f ). who performed a radiation magnetohydrody- 
namic (MHD) simulation of a stratified shearing box 
resulting in an average midplane radiation to gas pres- 
sure ratio of 14. Despite this large ratio, the simula- 
tion produced no obvious thermal instability over eight 
thermal times. The robustness of this result is some- 
what suspect, however, given that the computational 
methods employed precluded the photospheres (top and 
bottom) from being included within the simulation do- 
main. Moreover, 27% of the work done on the box disap- 
peared due to numerical losses, and half the mass was lost 
from the simulation domain during a heating/expansion 
phase. 

More recently, we ([Hirose. Krolik. fc Stonl [20060 de- 
veloped numerical techniques that conserve energy with 
high accuracy, place both photospheres within the prob- 
lem volume, and retain nearly all the initial mass on 
the grid. Solving the total energy equation enables us 
to capture all grid scale losses of magnetic and kinetic 
energy and convert them to internal energy in the gas. 
Furthermore, by simultaneously solving the internal en- 
ergy equation, we can compute the instantaneous local 
losses of magnetic and kinetic energy, i.e. the energy 
dissipation rate. The only violations of energy conser- 
vation come from the small amounts of energy (gener- 
ally ~ 0.1% of the total dissipated energy^) artificially 
injected by the action of the density fioor and similar 
numerical limiters. Using this code, we have previously 
studied a case in which gas pressure d ominated over ra- 
diation pressure ([Hirose. Krolik. fc Stone 2006.) and, in 
[Krolik. Hirose. fc BlaesI ([2007[ ): iBlacs. H irose. fc Krolik[ 
( 20071) . one in which the gas and radiation pressures were 
comparable. 

1 0.09% in "1112a" and 0.15% in "1126b"; see section[23]for the 
labels. 
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Both of the simulations studied in our previous pa- 
pers achieved a thermal balance between heating and 
cooling on time scales comparable to the thermal time. 
However, no real steady state was achieved in the simu- 
lation with comparable gas and radiation pressure. In- 
stead, there were long term (~ 7 thermal times) up and 
down variations in total energy content by factors of 
three to four. The two hottest epochs in this simula- 
ti on marginally violate d the t hermal instability criterion 
of IShakura fc SunvaevI (|1976D . but no thermal runaway 
occurred. 

In this paper we present results of two simulations 
that are in the fully radiation pressure dominated regime, 
with box-integrated radiation to gas pressure ratios ~ 10. 
The two simulations shared the same parameters, bound- 
ary conditions, and (almost) the same initial conditions, 
differing only in the small noise fluctuations imposed on 
the initial state to seed the MRI. Despite being in gross 
violation of the standard thermal instability criterion, 
there is no evidence for such an instability in either of 
their time evolutions, which in both cases extended over 
~ 40 thermal times. 

This paper is organized as follows. In section 2 we 
briefly summarize how the code works and discuss the 
parameters of the new simulations. Quantitative results 
about the thermodynamics, internal vertical structures, 
and variability of the simulations are presented in section 
3. In section 4 we discuss a simple model that might ex- 
plain the thermal stability we observe, and we summarize 
our findings in section 5. 

2. CALCULATION 

2.1. Basic Equations and Numerical methods 

The basic equations used in our simulations are the 
frequency-averaged equations of radiation MHD in the 
flux-limited diffusion (FLD) approximataion: 

dp 
dt 

d{pv) 



dt 



V • (pv) = 

V • (pvv) = -V(p + (?) 



(3) 



(V xB)xB + ^ '-^F + fsB 

Att c 



(4) 



de 

— +V-{ev)^~{S/-v){p + q) 



~{At:B — cE)k.qP — cEkcsP 



J diss 



— + V • (Ev) = -Vv : P + UttB - cE)Rlp 
at 

4fcB(r-rrad) „ „ 

+cEkosP ^ 5 - V • F 

dB 

— -Vx{vxB)^0 



cX 



(^ff ^cs)p 



VE 



(5) 

(6) 
(7) 
(8) 



The quantities p, e, and v are the density, internal 
energy, and velocity field of the gas. The pressure of 
the gas p is related to the internal energy by p — (7 — 
l)e, with adiabatic index 7 (here 7 = 5/3 is assumed). 



The quantity q is the pressure associated with a small 
artificial bulk viscosity. The quantity Qdiss represents 
the rate at which the internal energy must be changed 
fo r the total energy to be cons erved (see Appendix A 
in lHirose. Krolik. fc Ston^l2006[ ). The magnetic field is 
denoted by B. 

The radiation field is described by the energy density 
E, flux F, and pressure tensor P. In the FLD approxima- 
tion, the equations are closed with the relation P = fE. 
The Eddington tensor f is a function of the dimension- 
less opacity parameter R = \WE\/[{kg + Kcs)pE] and 
the flux limiter A. Here we adopt X{R) = (cothi? — 
1/R)/R (Levermore & Pomraning 1981). A more com- 
plete description of th e FLD approximation is found in 
iTiirner fc Ston^ (l20?)lh . 

The gas and radiation exchange both momentum 
and energy via electron scattering and free-free emis- 
sion/absorption. For free- free absorption, the Planck- 
mean opacity and Rosseland-mean opacity are com- 
puted as Kg = 3.7 X 1053(^9/g7)i/2 cmVg and = 
l.Ox 10^^(p^/e'')^/^ cm^/g. The electron scattering opac- 
ity is assumed to be Kos = 0.33 cm^/g. Energy ex- 
change via Compton scattering is represented by the 
third term in the RHS of the energy equations, where 
T and Tj-ad = {E/a Y are the gas and radia tion temper- 
atures, respectively (jBlaes fc Socratesll2003[ ). The quan- 
tity B in the free-free energy exchange terms is defined 
as caT^/ (47r). 

To simulate the dynamics in a section of an orbiting 
disk, we adopt the shearing box approximation, in which 
the sum of the Coriolis force, the tidal force, and the ver- 
tical gravitational force, /sb = —2p{^z) xv + Spil^xx — 
pH.'^zz, is added to the RHS of the momentum equation. 
Neglecting orbital curvature, we denote the radial coor- 
dinate, azimuthal coordinate, and vertical coordinate by 
the Cartesian coordinates x, y, and z. 

We solved these equations using a modified version 
of th e ZEUS code for radi ation MHD described in 
Hirose. Krolik. fc Sto^((200l . This version differs from 
the standard implementation of ZEUS in two main re- 
spects. First, strict energy conservation is achieved by 
capturing all numerical dissipation of kinetic and mag- 
netic energy, Qdiss, and delivering it to the gas's inter- 
nal energy. Second, the radiation diffusion equation is 
solved in time-implicit fashion using a multi-grid algo- 
rithm. Bot h of these alterations are desc ribed in Ap- 
pendix A of lHirose. Krolik. fc Stonl (|2006l ). 

2.2. Parameters 

Two parameters determine the problem in this context: 
the angular velocity VL and the surface density S. Using 
standard orbital mechanics and the radiation-dominated 
limit of the a model, these quantities can be written as 
functions of the central black hole mass M, the radius r, 
and the mass accretion rate M for a guessed value of a 
(Krolik 1999): 



GM 



3aKos M/Me \GM/6 



(9) 



(10) 
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As we have already explained, in a time-steady disk, the 
assumptions of radiative support against vertical gravity 
and energy conservation can be used to relate the disk 
vertical thickness H and the radiated flux Fq to M and 



Fn = —MVL'- 



H = 



3 KesM 

Stt c 



(11) 



(12) 



However, properly speaking, H and i^o (as well as M) are 
driven by local disk dynamics and are in general time- 
dependent; thus, these relations should be regarded as 
guesses that may or may not be confirmed. For this 
reason, we use the a model prediction for H merely as 
an order-of-magnitude estimator of the actual vertical 
thickness, and adopt it as our unit of length. 

The specific parameters we chose for the simulations 
reported here are summarized in Table 1. The efficiency 
•q is assumed to be 0.1. 

2.3. Grid and Boundary Conditions 

The size of the computational domain is {L^, Ly, Lz) — 
{0A5H,1.8H,8AH). It is divided into (48,96,896) cells 
with a constant cell size Ace = Az = Ay/2. The az- 
imuthal boundaries are periodic and the radial bound- 
aries are shea ring-periodic in the local shcaring-box app- 
proximation (|Hawlev. GammieTfc Balbus 1995). The 
vertical bou ndaries are (with one exc e ption ) treated ex- 
actly as in iHirose, Krohk. fc Stond (|2006f) : hydrody- 
namic quantities are given outflow (free) boundary con- 
ditions, while the diffusive radiation flux and the elec- 
tric field do not change across the boundary. We refer 
the interested reader to the lengthy discussion found in 
that paper regarding the difhculties of treating a sub- 
sonic boundary and the requirement for introducing a 
small resistivity in the boundary cells. The only dif- 
ference between our treatment and the one described 
there is that the artificial resistivity is extended into 
the problem volume, tapering smoothly from the level 
in the ghost cells to zero at 3 2 cells from the boundary 
(jKrolik. Hirose. fc Blaesl[200l . 

2.4. Initial Conditions 

The disk segment in the simulation box is initially as- 
sumed to be in both dynamical and thermal equilibrium 
in the vertical direction, ignoring any magnetic forces. 
The key assumption underlying the structure of the ini- 
tial cond ition is the verti cal profile of the dissipation rate; 
foUowing'Hiro se. Krolik. fc Stonj (|2006f ). we assume that 
Qdiss is proportional to p/\/t, where r is the optical 
depth from a given point to the nearest surface. The 
following equations are then solved to compute the ver- 
tical profiles of gas and radiation field for 1 ^ t ^ tq, 
where tq is the optical depth to the midplane. It is re- 
lated to the surface density S by tq = (S/2)Kes + 1; for 
our parameters, tq = 8.79 x 10"^. 



dp 

dz 
d7 



_P_ 

RT 



ZpF 



1 



(13) 
(14) 



dE 

dF_ 
d7 



3F 



diss 



^csP 



(15) 
(16) 



Our procedure is to begin with the fourth equation, for 
which we impose the boundary conditions F(to) = 
and F(l) = Fq; the result of solving this equation is 
the flux profile F{t). Then the radiation energy profile 
E{t) is obtained immediately from the third equation. 
Next, the first two equations are solved numerically with 
boundary conditions z{to) = and p{l) = pfloor, to get 
the density profile p(r) and the (inverse) opacity profile 
z(t). The density floor pfloor is set to 10~^p(ro). The gas 
temperature T is assumed to be equal to the radiation 
temperature [Eja)^!^. Then the midplane density and 
the height of the disk surface are obtained from p(ro) = 
5.66 X 10-2 g cm-2 and z(l)/i? = 1.36. Outside the 
photosphere [z > z(l)), we assume a uniform gas density 
Pfloor, a gas temperature r(z(l)), and a radiation energy 
density £^(^(1)). 

The initial velocity field is determined so that the in- 
ertial force is zero {—2p{ilz) x v + Spil^xx = 0), i.e., 
it is only the orbital shear v = (0, —3nx/2, 0). The ini- 
tial configuration of the magnetic field is a twisted az- 
imuthal flux tube (with net azimuthal flux) placed at 
the center of the simulation box. The field strength in 
the tube is uniform and the ratio of the poloidal field to 
the total field is 0.25 at maximum. The flux tube cross 
section is an oval with radius in x of 0.75(La;/2) and ra- 
dius in z of 0.80z(l). The initial box-integrated plasma 
P = {{E/3) +p)/iB'^/8Tr) = 11.5, while in the midplane 
in the initial state it is 40. The corresponding maximum 
MRI wavelength in the midplane is resolved with 18.7 
grid zones. 

The calculation is begun with a small random pertur- 
bation in the poloidal velocity. The maximum amplitude 
of each velocity component is 1% of the local sound speed 
defined as ^y {{A / 3) {E / 3) + ■yp)/p- Two nearly identical 
simulations were run for just over 600 orbits each, the 
only difference between them being the seed for the ran- 
dom initial poloidal velocity perturbations. We refer to 
these two simulations as 1112a and 1126b, and we now 
discuss their properties. 

3. SIMULATION RESULTS 
3.1. Technical Consistency 

We are interested in the stability of disk segments in 
which the surface mass density is constant. To keep E 
steady, the box height must be chosen carefully because 
our boundary conditions permit outflows at the top and 
bottom of the simulation box, but forbid inflows. Con- 
sequently, mass is lost comparatively readily in a shorter 
box, but mass can be added in a taller box due to more 
frequent application of the code's density floor when in- 
flows are shut off. In simulations 1112a and 1126b, we 
set the box height at 8AH, which results in variations of 
the surface mass density of less than 2.5% (FiglJ). 

We must also ensure that the MRI is well-resolved 
throughout the box at all times. The fastest growing 
mode has a vertical wavevector and wavelength A* = 
6A9/il^y B^/ p. Averaging A* over horizontal planes, we 
show the number of grid-cells per wavelength of this 
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TABLE 1 



Parameters 


Definition 


Vaiuc 


Comment 


M 


6.62 Mq 


1.31 X W'-'^g 


Mass of central black tiole 


r 


30 (GM/c2) 


2.93 X let'em 


Radius 


M 


0.1{Le/c^)/v 


7.52 X lO^^'g 


Accretion rate 


a 


0.0125 


Stress/pressure ratio 


n 


eq.m 


1.90 X 10^ 


Orbital frequency 


Fo 


eq.El] 


4.62 X 10^1 erg cm-^ s'l 


Surface radiation flux 


E 


eq.[TO] 


5.32 X 10'' g cm-2 


Surface density 


H 


eq.[ll 


1.46 X 10^ cm 


Disk scale height 





4 - 




-2 - 



-4 - 



100 200 300 400 500 600 

t (orbits) 

Fig. 1. — Variation (in percent) of surface mass density as a 
function of time (S(t) - S(i = 0))/E{i = 0) in 1112a(black) and 
1126b{gray). 



100 2«) 30O .100 300 600 

100 200 30O 4O0 500 600 



Fig. 2. — Top panel: Number of cells per wave length of the most 
unstable mode of MRI in a horizontally averaged sense (color) and 
time variation of the photosphere location defined as the Eddington 
factor / = 0.5 (solid curves) in 1126b. Bottom panel: Same as the 
top panel except for 1112a. 



mode in Fig. [5] There are at least 8 cells per wavelength 
(and generally a great many more) everywhere in the box 
except near the midplane during a few brief epochs (e.g. 
r~ 400 in 1112a). 

The quality of our results also depends on keeping both 
the top and bottom photospheres within the simulation 
domain at all times. Fig. [2] also shows the locations 
of the photospheres in both simulations; we define the 
photosphere by the place where the Eddington factor 
f{R) = 0.5, for R the horizontally-averaged dimension- 
less opacity parameter. The photospheres are comfort- 
ably within the computational domain throughout the 
duration of both simulations. 



3.2. Thermal Balance 

The thermal history of the two simulations is shown in 
Figs. [3] and [H The top panel in each figure shows the 
volume integrated dissipation rate per unit surface area 
of the box (red) and the total radiative flux leaving the 
top and bottom faces of the box (blue), both on a lin- 
ear scale. This panel shows clearly that the box is able 
to maintain a close thermal balance on timescales very 
short compared to the durations of these simulations. 
The bottom panels show the total radiation, gas inter- 
nal, magnetic, and turbulent kinetic energies in the box 
on a logarithmic scale to emphasize better the relative 
fluctuations in each of these quantities. If we define the 
unit of energy content as the gas thermal energy and con- 
sider both simulations together, the radiation energy can 
range from ~ 5-30, the magnetic energy from ~ 0.15- 
1, and the kinetic energy from ~ 0.05-0.25. Adopting a 
nominal radiative efficiency of 0.1, we find that the mean 
stress in simulation 1112a would drive an accretion rate 
relative to Eddington of 0.17, while in 1126b it would 
have been 0.23. Thus, even when averaging over rela- 
tively long timescales, the mean accretion rate for this 
surface density and orbital frequency is still uncertain at 
the level of several tens of percent. 

A number of interesting results are apparent in these 
figures, but the one that is most significant is the fact 
that an approximate thermal equilibrium has been es- 
tablished in both simulations. The thermal time as mea- 
sured by the radiation and gas internal energies divided 
by the total radiative flux has short time scale fluctua- 
tions, but averaging over timescales ~ 50 orbits, we find 
that in both simulations the instantaneous cooling time 
can be anywhere from ~ 10 orbits at times of relatively 
low energy content to ~ 20 orbits at times of especially 
high energy in the box. Averaged over the entire dura- 
tion, the cooling time is 13 orbits in 1112a, and slightly 
longer, 15.5 orbits, in 1126b. Thus, both simulations ran 
for ~ 40 thermal times, and yet there is no evidence 
of run away heating and cooling. IShakura fc SunvaevI 
(|1976l ) predicted that thermal instability should occur 
whenever the ratio of radiation to gas pressure exceeds 
3/2, but the radiation to gas energy ratios encountered 
here translate to pressure ratios ~ 2.5-15, well above the 
predicted instability threshold. 

Although it does appear that a crude long-term mean 
value for the energy content can be estimated, the ampli- 
tude of fluctuations on timescales a significant fraction of 
the simulation duration is still quite large. We quantify 
these fluctuations by several measures. 

First, we note that transients associated with the ini- 
tial condition are generally erased by ~ 100 orbits into 
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3 1020 
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^^^^^^ 

100 200 300 400 500 600 

t (orbits) 

Fig. 3. — Top panel: Total volume integrated dissipation (red) 
and radiative cooling (blue) rates as a function of time in simulation 
1112a, plotted on a linear vertical scale. Bottom panel: Total 
energy content in the box as a function of time in simulation 1112a, 
in the form of radiation energy (red) , gas internal energy (green) , 
magnetic energy (blue), and turbulent kinetic energy (black). (The 
turbulent kinetic energy is defined as the kinetic energy associated 
with the fluid velocity measured relative to the background shear 
flow.) Note that the bottom panel is plotted on a logarithmic 
vertical scale. 
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Fig. 5. — Variability power spectral densities of volume- 
integrated magnetic energy (black) and radiation energy (gray) in 
simulation 1112a. Both time-dependent energies were normalized 
by their time-averaged values before computing the power spectral 
densities, and both power spectra are plotted after having been 
multiplied by the frequency. The vertical dashed curve indicates 
the expected frequency of the longest wavelength vertical acous- 
tic wave. Dot-dashed curves are broken power law fits to the two 
power spectra in the frequency range 0.01 to 10.0, scaled by inverse 
orbital periods. 
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Fig. 4. — Same as Fig. \3\ except for simulation 1126b. 



the simulation. Discounting that initial period, we see 
that the typical fluctuation amplitudes in the radiation, 
magnetic, and kinetic energy are similar, all factors of 
several, but the fluctuation levels in the gas energy are 
much smaller, only 6% rms. This means that the dissi- 
pated magnetic and kinetic energies are readily converted 
to radiation energy. On the other hand, although there is 
considerable high-frequency power in the magnetic and 
kinetic energy fluctuations, there is much less in the ra- 
diation energy. That this should be so is not surprising: 
the turbulence driving the magnetic and kinetic energy 
densities has a characteristic timescale of order an or- 
bit, whereas the radiation energy changes on a thermal 
timescale, an order of magnitude longer. 

These remarks can be quantified by studying the power 
spectra of the fiuctuations (Fig. [5] illustrates the case of 
1112a; 1126b is qualitatively similar). The offset between 
the magnetic and radiation Fourier power spectra is a 
measure of the larger contribution from short timescale 
fluctuations in the time series of magnetic energy. Only 
at the very longest timescales does the variance in the ra- 
diation energy become similar to that in the magnetic en- 
ergy (and also at very short timescales, where the fluctua- 



tion amplitude is extremely small). More quantitatively, 
both power spectra can be approximately described by 
broken power-laws: 



Pif) 



^-2.38, 



-3.65. 



/break — 0.171 
-^•'^ /break = 0.118 



magnetic 
radiation ' ^ ' 



The first number in each exponent is the low-frequency 
slope; the second number is the high-frequency slope. 
In both cases, the break occurs at ~ 5-10 orbital pe- 
riods, about half the thermal timescale. At high fre- 
quencies, both power spectra decline quite steeply. On 
the other hand, at low frequencies, both power spec- 
tra are "red" enough that the variance is dominated by 
the low-frequency cut-off posed by the duration of the 
simulation. This "infrared divergence" is weak — almost 
logarithmic — for the magnetic fluctuations, but is rather 
stronger for the radiation spectrum. As we will argue in 
§ 4, the significant power at very long timescales in the 
magnetic fluctuations is an important element in explain- 
ing these disk segments' thermal stability. As we will 
also argue, the fact that the two power spectra approach 
one another at the lowest frequencies is a symptom of 
the fact that MHD turbulence drives both, so that they 
must approximately coincide on timescales longer than 
the thermal equilibration time. 

The essence of the a model is the expectation that 
the stress should be proportional to pressure. We exam- 
ine the validity of this expectation in Fig. [6l but using 
the box-integrated magnetic energy as a stand-in for the 
stress. As this figure shows, the two do tend to vary to- 
gether, but with order unity fluctuations about the trend. 
The logarithmic slope of the least squares fit shown in 
the figure is 0.71, indicating that the relation is shal- 
lower than linear. For future reference (this relationship 
will be called upon in § 4), we also display (Fig. [7|) the 
analogous relationship between tcooi and the integrated 
radiation energy. The logarithmic slope of the fit in that 
case is 0.32. It is worth pointing out that the slopes of 
these relations are subject to some chaotic variation: the 
magnetic-radiation energy correlation has slope 0.79 and 
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Fig. 6. — Box-integrated magnetic energy and radiation energy 
compared every 0.01 orbit from T = 100 to T = 600 in simulation 
1112a. Both energies are in units of lO^'^ erg cm~^. The dashed 
hne is a least squares fit in the logarithms of these two quantities. 



1.0 



0.8 - 




0.2 - 



0.0 1 

-40 -20 20 40 

lag (orbits) 

Fig. 8. — Cross correlation with volume integrated magnetic 
energy of volume integrated radiation energy (solid), gas internal 
energy (dotted), and turbulent kinetic energy (dashed) as a func- 
tion of time lag for simulation 1112a. Negative lag implies that the 
energy in question lags behind the magnetic energy. 




Radiation Energy 

Fig. 7. — Thermal time (in units of orbits) and radiation energy 
(in units of lO^'' erg cm""^) compared every 0.01 orbit from T = 100 
to T = 600 in simulation 1112a. The dashed line is a least squares 
fit in the logarithms of these two quantities. 

the thermal time-radiation energy correlation has slope 
0.44 in simulation 1126b. 

Casual study of Figs. [Sj HI and [6] suggests that, after 
allowance for the smoother variations in the radiation 
energy, all the contributions to the energy content are 
well-correlated in time. This is almost true. A cross- 
correlation between them (Fig. [5]) reveals, however, that 
fluctuations in the magnetic energy lead those in the radi- 
ation energy (and those in the gas energy) by 5-15 orbits, 
roughly a thermal time. In § 4, we will show that this 
fact is a crucial clue for understanding the boxes' ther- 
mal stability. The turbulent kinetic energy, on the other 
hand (and unsurprisingly) is closely coincident in time 
with the turbulent magnetic energy. 

We close this subsection by remarking on the quasi- 
periodic oscillation (QPO) that appears in the power 
spectrum of radiation energy at a frequency slightly 
higher than the orbital frequency. We believe this to be 
a vertical acoustic oscillation. Standing vertical acous- 
tic waves in a polytropic, Newtonian, Keplerian disk can 
occur at angular frequencies given by 



n:,{n^ + 2n - 1) 
2n 



1/2 



(18) 



where is an integer greater than or eq ual to 

two and n is the polytropic index (|Katol 120051 : 



iBlaes. Arras, fc Fragile! [20061 ) . The lowest frequency of 
this spectrum, nz = 2, fits the QPO exactly for a ra- 
diation pressure supported configuration (n = 3). The 
eigenfunction of this mode takes the form of a vertical 
breathing mode, with fluid displacements expanding and 
contracting symmetrically about a node at the midplane. 
As we discuss in section 3.4 below, this mode plays an 
important role in reconciling the actual vertical profiles 
of turbulent stress and dissipation with the constraints 
of thermal and hydrostatic equilibrium. 

3.3. Stress 

In the original form of the a-model, it is supposed that 
the stress bears a fixed ratio to the total pressure. How- 
ever, it has also b een proposed (ISaki moto fc Coronitil 
119811: iTaam fc LiiJ [19841. usually with a view toward 
quenching the thermal instability, that the stress might 
instead scale more closely with the gas pressure alone, or 
with the geometric mean of the gas and radiation pres- 
sures. We can test these hypotheses with our simulation 
data. 

Fig. [5] shows the time dependence of the box aver- 
age of the r0 component of the Maxwell and turbulent 
Reynolds stress in simulation 1112a. The fact that it 
has considerable short timescale fluctuation power im- 
mediately rules out the possibility that it should respond 
to either the gas pressure or the radiation pressure in 
all respects because their histories show only very weak 
short timescale fluctuations. The large amplitude fluc- 
tuations in the stress over long timescales similarly show 
that it cannot be driven, even in an averaged sense, by 
the gas pressure, for the fluctuations in the gas pressure 
are far too small. That the volume-averaged stress is at 
all times well below the critical stress for pure radiation 
support demonstrates that, although radiation pressure 
truly dominates in the central part of the disk, other 
mechanisms (notably magnetic pressure, as we will show 
in § 3.4) account for support against vertical gravity in 
substantial portions of the total simulation volume. 

To make this point very explicit, we show in Fig. (TU] 
the time-dependence of the ratio of the box-integrated 
stress to the different box-integrated candidate pressure 
quantities for simulation 1112a. That is, these curves 
show the time-dependent Shakura-Sunyaev a parameter 
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Fig. 9. — Time dependence of the box-averaged r(j> stress in 
simulation 1112a. The horizontal dashed line indicates the 2cQ/3k 
stress expected from hydrostatic and radiative equilibrium. 
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Fig. 10. — Time dependence of the ratio a of the box-averaged 
r<f> stress to various box-averaged pressures in simulation 1112a. 
From top to bottom, the stress is scaled with gas pressure, the 
geometric mean of gas and radiation pressures, and with total (gas 
plus radiation) pressure, respectively. 

for various stress prescriptions. The best would presum- 
ably be the one that gives a value of a that varies least 
with time. The least variation (taking the stress divided 
by the gas plus radiation pressure) is a factor of 4, while 
the greatest variation (found when the stress is measured 
in units of the gas pressure) is a factor of 10. In no case, 
therefore, does the stress follow exactly the fluctuations 
in any of these pressure measures. 

It is worth remarking in this context that the mean 
value of the ratio of integrated stress to integrated total 
pressure was 0.023 in 1112a and 0.019 in 1126b. In our 
previous simulations, this number was 0:^ 0.03 (for the 
Pr Po simulation reported in Kr olik. Hirose. fc BlaesI 
(120071)) and 0.016 (for t he pj pn ~ 0.2 simulation of 
iHirose. Krolik. fc Stond (|2006f )). Thus, we see little 
trend in this ratio as a function of the ratio of radia- 
tion to gas pressure. In addition, this series of simu- 
lations is, somewhat inadvertently, a crude test of nu- 
merical convergence. Due both to the requirements of 
the simulations and to the increase in available com- 
putational power over time, the resolutions employed 
have improved steadily from one to the next. Measured 
in terms of our estimated scale-height, the cell-size Az 
changed from 0.06257f in the gas-dominated simulation 
to 0.0234iJ in the pr ~ Pg simulation to 0.0094iJ in 
the simulations presented here. The rough constancy of 



the stress/pressure ratio across these simulations is of 
some interest in view of reports that unstratified shear- 
ing box simulations show a dec hning ratio of stress to 
press ure as resolution improves (jFromang &: Papaloizod 
120071. It should be pointed out, however, that these un- 
stratified simulations had zero net magnetic fiux, whereas 
our stratified simulations have a net azimuthal fiux, al- 
beit one that is not preserved by the boundary conditions 
and fluctuates in direction over the course of the simula- 
tions. 

3.4. Vertical Energy Transport and Hydrostatic 
Balance 

As first pointed out bv lShakura fc SunvaevI (|1976l ). in 
any radiation-dominated geometrically-thin disk, there is 
a characteristic emissivity per unit volume: = cfl^ /k. 
Its existence follows immediately from attributing hydro- 
static balance to radiation force. When that condition is 
met (and the opacity is purely electron scattering), 

K,,F,/c = zfl^; (19) 

differentiating both sides with respect to height z gives 

^^ = ^j = f^2 (20) 
c dz c 

If all the dissipation is transformed locally into photon 
energy, there is a corresponding characteristic dissipation 
rate = j*. We find (Fig. [TT]) that the dissipation rate 
is comparable to through most of the disk body, but 
is typically somewhat greater at \z\ ~ 0.5H (by ~ 30%) 
and rather less in the upper disk. At all altitudes it is 
dominated by magnetic losses. 

If hydrostatic balance obtains (as we will show shortly, 
departures are almost always small), where Q > Q* over 
a significant range in height, it must be that j < Q be- 
cause otherwise the fiux would exceed the value whose 
force balances gravity. Consequently, not all the energy 
liberated by dissipation goes directly into radiation fiux. 
As shown in Fig. [T^l for \z\ < 2H, a noticeable (but 
minority: ~ 10%) fraction of the vertical energy fiux is 
in the form of photons adverted upward with the mat- 
ter (Evz), rather than diffusing through it (Fj); because 
this component does not move through the matter, it 
exerts no force. Like the total radiation content in the 
box, the advective energy fiux also oscillates with the fre- 
quency of the fundamental vertical breathing mode. If 
this oscillation behaved with exact symmetry in its out- 
ward and inward moving phases, it would carry no net 
energy. However, at the outer surface of the region of 
significant advection {\z\ ~ 27?), we see a divergence of 
the diffusive flux that exceeds the local dissipation rate 
(Fig. [in]); that is, radiation energy diffuses out of the up- 
welling matter at the top of its oscillatory range. In this 
respect, the radiative advection seen here resembles con- 
vection, although the vertical gas motions are acoustic in 
nature and not due to buoyancy. In addition, a smaller 
amount of energy is carried upward by work done by ra- 
diation pressure forces and by Poynting flux {E x B)^, 
which also oscillates with the frequency of the breathing 
mode. Clearly, MHD flux-freezing causes magnetic en- 
ergy to be advected with the gas along with radiation 
energy in the breathing mode. The energy transported 
by radiation advection, work done by radiation pressure, 
and Poynting fiux is deposited (or dissipated) around 
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Fig. 11. — Horizontally and time-averaged total dissipation rate 
Q (black) and stress times rate of strain (dotted) in simulation 
1112a. The total dissipation rate is the sum of magnetic energy dis- 
sipation (red) and kinetic energy dissipation (blue). The horizontal 
dashed line marks the characteristic dissipation rate Q, = cQ'^ /k. 

H ^ \z\ < 2H. Above this point, the energy flux is over- 
whelmingly due to radiation diffusion. Energy transport 
by gas advection (e -I- pv^/2)vz and gas pressure work is 
negligible throughout the simulation domain. 

Fig. [11] also shows the vertical profile of stress times 
rate of strain, — Tr^dfJ/dlnr. This profile is significantly 
more concentrated near the midplane and at z = ±0.5H 
than the dissipation profile, indicating that the work 
done by the turbulent stresses is, in the time average, 
partly dissipated locally and partly converted into me- 
chanical energy of gas motion, which is then transported 
vertically primarily by Poynting flux and radiation pres- 
sure work associated with the breathing mode and dissi- 
pated farther out. This energy is deposited well outside 
the peak region of stress, and subsequently carried fur- 
ther outward by radiative diffusion. This supplemental 
energy transport mechanism contributes to dFz/dz, help- 
ing to maintain radiation-aided support against gravity 
at altitudes higher than one might expect if dissipation 
correlated perfectly with local stress. 

The net result of these processes is illustrated most 
clearly in Fig.[T31 As this figure shows, despite the peaks 
in Q that reach well above Q*, the divergence of the ra- 
diation flux Fz is kept just under across the entire 
central body of the disk, from z ~ —1.5/7 to z ~ +1.5H. 
It is the radiation advection flux Evz that compensates 
for the excess dissipated energy (Q — Q*) and carries it 
away. Thus, thermal balance is locally maintained in the 
entire region, as evidenced by the fact that the divergence 
of the summed radiation and gas energy fluxes (Fz+Evz) 
matches the dissipation rate Q very well. At the same 
time, Poynting flux {E x B)z and the flux of work done 
by radiation pressure forces transports non-dissipated ex- 
cess mechanical energy away from the midplane to be 
dissipated further out. The sum of the divergences of 
all these energy fluxes (together with advection of gas 
energy and work done by gas pressure, which are both 
small) compensates exactly for the work done on the fluid 
by turbulent stresses, thereby establishing global energy 
conservation. 

Consistent with the rough match between Q and Q*, 
most of the vertical support, particularly in the center 
of the disk, is due to radiation forces (Fig. [Ti|) and hy- 
drostatic balance is maintained quite closely. In fact, 
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Fig. 12. — Horizontally and time-averaged diffusive radiation en- 
ergy flux Fz (solid red), advected radiation energy flux Evz (dashed 
red), Poynting flux (E X B)z (blue), and advected gas energy flux 
{e -\- pv'^ /2)vz (green) in simulation 1112a. 
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Fig. 13. — Horizontally and time-averaged vertical derivatives 
of diffusive radiation energy flux Fz (solid red), advected radiation 
energy flux Evz (dashed red), Poynting flux {E X B)z (blue), and 
advected gas energy flux {e + pv^ /2)vz (green) in simulation 1112a. 
In addition, the grey curve shows the profile of P : — xpF ■ v/c, 
which is equivalent to the divergence of the flux Evz/3 of work done 
by radiation pressure when the medium is optically thick. The sum 
of the diffusive and advective radiation energy flux derivatives is 
shown as the solid black line, and matches well the horizon tally 
and time-averaged dissipation profile Q (solid black line in Fig. llll l. 
The sum of all the flux derivative profiles is shown as the dotted 
black line, and matches well the horizontally and time-aver aged 
profile of stress times rate of strain (dotted black line in Fig. [TTJ. 
For reference, the horizontal dashed line marks the characteristic 
dissipation rate Q* = cH^/k. 

departures from time-averaged hydrostatic balance are 
smaller, in fractional terms, than any local deviation 
(again time-averaged) between Q and Q^. In general 
terms, this occurs because the timescale for dynamical 
equilibration (~ 1/S^) is the shortest timescale in the sys- 
tem. In detail, two different effects account for the fact 
that deviations from hydrostatic balance are so small. As 
just discussed, in the disk body, where Q > Q<,, some of 
the vertical energy transport is carried by other mecha- 
nisms, notably radiation advection (Fig. [T2)) . At higher 
altitudes, vertical support depends more on the gradient 
of magnetic pressure than on the gradient of radiation 
pressure. 

The structure that is established by these several 
forces is very similar qualitatively to the structure seen 
in all previous vertically-stratified shearing-box simula- 
tions with genuine thermodynamics: a central region 
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Fig. 14. — Horizontally averaged vertical accelerations at the 
200 orbits epoch in simulation 1112a. The different curves show 
accelerations from radiation pressure (red), gas pressure (green), 
and magnetic forces (blue). The total of these accelerations is the 
black curve, and is comparable to the straight black line indicat- 
ing the gravitational acceleration. The dashed blue curve is the 
contribution to the magnetic acceleration from magnetic pressure 
gradients, while the dot-dashed blue curve indicates the contribu- 
tion from magnetic tension. The vertical dashed and dot-dashed 
lines indicate the locations of the scattering and thermalization 
photospheres of the horizontally-averaged structure, respectively. 
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Fig. 15. — Horizontally and time-averaged vertical density profile 
in the 1112a (solid) and 1126b (dashed) simulations. 

{\z\ < 3H) supported by a combination of gas and ra- 
diation pressure lying within an outer region dominated 
by magnetic forces. Within the magnetically-dominated 
region, the outward force due to the magnetic pressure 
gradient is substantially cancelled (in the mean) by the 
inward force due to magnetic tension (Fig. [T4|) . 

The resulting profiles of density and pressure are 
roughly exponential (Figs. [15] and I16|) . An exponen- 
tial description is particularly good for the density, with 
time-averaged scale-height ~ 0.8H. However, we cau- 
tion that the instantaneous scale-height can vary by or- 
der unity as the radiation energy varies up and down. It 
is also appropriate to remark at this point that, contrary 
to the simple assumption made in classical a-models, the 
density is very far from constant because the dissipa- 
tion per unit mass is also very far from constant. As 
we have found in earlier vertically-stratified studies with 
smaller ratios of radiation to gas pressure, the dissipa- 
tion rate per unit mass increases with decreasing mass 
column density to the nearest outer surface. 

The pressures, while still behaving roughly as expo- 
nentials in height, have additional structure. The radia- 
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Fig. 16. — Horizontally and time-averaged vertical profiles of ra- 
diation pressure (red), gas pressure (green) and magnetic pressure 
(blue) in the 1112a (solid) and 1126b (dashed) simulations. 
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Fig. 17. — Horizontally and time-averaged gas temperature 
(solid) and radiation temperature (dotted) in 1112a (black) and 
1126b (gray). 

tion pressure declines more steeply within \z\ < 3H than 
without because the density drops so steeply outward. 
From the photosphere outward, the radiation pressure 
changes only slowly with height because that is the free- 
streaming limit in plane-parallel geometry. The magnetic 
pressure profile, on the other hand, is roughly flat-topped 
in the central disk, with a pair of weak local maxima at 
\z\ ~ H. In that central region, the plasma f3 ~ 50, but it 
drops rapidly outward, falling below unity at \z\ ~ 2.5H 
where the magnetic pressure begins to dominate the to- 
tal. 

Fig. [T7] shows horizontally and time-averaged gas tem- 
perature T and radiation temperature T^ad = (E/aY^^, 
while Fig. [TSl shows the height of the scattering and ther- 
malization photospheres as functions of time for sim- 
ulation 1112a. These photospheres are defined as the 
surfaces where the horizontally averaged scattering op- 
tical depth and geometric mean of the scattering and 
Planck mean absorption optical depths, respectively, 
equal unity. The time-averaged height of the thermaliza- 
tion photosphere in simulation 1112a was 2.0H, in ap- 
proximate agreement with the height at which the time- 
averaged gas and radiation temperatures begin to sepa- 
rate in Fig. [T71 The scattering photosphere is of course 
further out, with a time-averaged height of 3.3iJ. 

Fig. [m shows that very large density fluctuations are 
always present at both the thermalization and scatter- 
ing photospheres of the horizontally averaged structure. 
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Fig. 18. — Horizontally-averaged height of the scattering (black) 
and thermalization (gray) photospheres as functions of time in 
1112a. 
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Fig. 19. — Histogram of the ratio of density to horizontally- 
averaged density at the thermalization (green) and scattering (red) 
photospheres of the horizontally averaged structure at epoch 200 in 
simulation 1112a. The upper and lower photospheres are indicated 
by solid and dashed histograms, respectively. The fluid motions 
are highly compressible in this region due to the large ambient 
magnetic forces. 

Similarly large density fluctuations in the photospheres 
have been seen in all our previous simulations (see Fig. 
16 in lBlaes. Hirose. &: Krolikll2007f ). and are due to the 
fluctuating magnetic forces that dominate all other forces 
in these regions. 

We close this section by remarking that we have con- 
firmed many of the results s een in t he ra diation pres- 
sure dominated simulation of [Turned (|2004D . despite the 
fact that that simulation was not fully energy conserv- 
ing, nor did it include the photospheres within the sim- 
ulation domain. In particular, he also found no thermal 
instability and a vertical density profile that was highly 
concentrated toward the midplane. As shown in Fig. I20| 
we also find similarly small density fluctuations near the 
midplane, of at most a factor ^ 2 between maximum 
and minimum. The magnetic pressure proflle in his sim- 
ulation had a double-peaked structure that was lower 
than the gas pressure at the midplane, but exceeded the 
gas pressure further out, also in agreement with what 
we have found. However, his magnetic pressure nowhere 
exceeded the radiation pressure, in contrast to the outer- 
most layers in our simulation. Moreover, he found that 
radiative diffusion carried only two-thirds of the outward 
energy flux, with radiation advection being the most im- 
portant secondary coolant. These differences can per- 
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Fig. 20. — Histogram of the ratio of density to horizontally- 
averaged density at the midplane for simulation 1112a at epoch 
200. The fluid motions are highly incompressible in this region. 

haps be ascribed to the fact that his simulation did not 
include the photospheres. If we neglected the data from 
our simulations outside \z\ = 2H, we would see results 
qualitatively similar to his. 

The o nly re spect in which we really disagree is that 
iTurneij (I2004D also claimed that radiation damping 
(jAgol fc Kroii5ll998[ ) contributed 29% of the heating in 
the midplane regions of his simulation. However, we have 
demonstrated a complete accounting of the thermal en- 
ergy balance without calculating the contribution from 
this process at all. Turner computed the radiation damp- 
ing rate by integrating — PV • v over the midplane region 
of his box, with P = E/3 where radiation dominates and 
is isotropic. This i s a reasonable met hod in a nonstrat- 
ifled shearing box ([Turner. Stone &: Sa no 2002.) because 
the volume is fixed: consequently, the net compressive 
work indicates non-adiabatic behavior. However, in a 
stratified shearing box, PV • v is really a piece of the di- 
vergence of the flux of work done by radiation pressure, 
V • (Pv), and we have shown that this plays a contribut- 
ing role to transporting excess mechanical energy out of 
the midplane by the breathing mode to be dissipated 
further out. Moreover, even if we follow Turner in using 
— PV • V to provide an estimate of the radiation damp- 
ing near the midplane, we then find that this can be at 
most 1.3% and 0.7% percent of the dissipation in the 
midplane region in 1112a and 1126b, respectively. It is 
conceivable that one of two limitations in Turner's code 
may be responsible for this quantitative difference be- 
tween his simulation and ours: the lack of a photosphere 
and the inability to capture grid scale losses of kinetic 
energy. Without performing additional simulations that 
separately reproduce each of these possibilities, we can- 
not be certain as to the exact origin of the discrepancy. 
Nevertheless, we suspect that lack of true energy conser- 
vation in Turner's simulation is the primary cause. 

4. INTERPRETATION 

In the traditional a-model approach, whose logical un- 
derpinning is dimensional analysis, the stress, and there- 
fore the magnetic field energy density, is scaled to the 
pressure. Considered in a logically rigorous fashion, di- 
mensional analysis suggests that two quantities with the 
same units are comparable in magnitude, but does not 
determine a causal relation between the two: that is, 
it does not tell us which quantity determines the other. 
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Nonetheless, if it is easier to estimate one than the other, 
it is easy to slip into the thought that the one harder to 
estimate follows the other, not just in scale, but in time. 
In this case, the a-model has often been imagined to 
suggest that the pressure controls the magnetic field. 

However, a crosscorrelation of the magnetic field en- 
ergy against the radiation energy in this simulation shows 
that fluctuations in the magnetic field lead those in the 
radiation. If so, the dominant sense of causality must 
be the other way around — the magnetic field controls 
the radiation, not vice versa. In fact, that is the direc- 
tion suggested by the dynamics: The MRI causes field 
fluctuations to grow, deriving the needed energy from 
orbital shear; nonlinear mode-mode interactions move 
power from long lengthscales to short; short wavelength 
dissipative effects transfer energy from the magnetic field 
to heat in the plasma; finally, electrons in the plasma 
radiate photons. In other words, intrinsic fluctuations 
in the magnetic turbulence can create a pressure-stress 
correlation even when there is no direct influence of the 
pressure on the saturation level of the magnetic field. 
All that is required is that their amplitude on timescales 
long compared to a thermal equilibration time are large 
enough, and these are already known to occur even in 
gas-dominated cas es in whi ch the pressure varies very 
httle (Hirose, Kro lik. fc Sto ne 2006). When that is the 
case, the increased dissipation rate due to the stronger 
field leads to higher pressure, but not necessarily to 
any subsequent increase in the stress. Thus, a key as- 
sumption of the argument for thermal instability in the 
radiation-dominated regime is undermined. 

Using an approach reminiscent of t he long - wavel ength 
limit of the formalism developed in iPiranI ()1978) . the 
qualitative picture outlined above can be modelled by 
the following pair of equations: 



E 



BO 



Ef 



dt ^growth \Ero 

dEji Eb Er 



- ^ (21) 
(22) 



dt ^diss ^cool 

where Eb^r are the box-integrated magnetic and radi- 
ation energies, Ebo,ro are their equilibrium values, and 
R{t) is a random function with expectation value unity 
describing the intrinsic fluctuations in the scaling of mag- 
netic energy density with pressure. We also suppose that 
the equilibrium magnetic energy density grows with pres- 
sure (here the radiation energy) with logarithmic deriva- 
tive n. Thus, we allow for feedback in the sense that the 
pressure might be able to influence the stress. Lastly, 
^growth is the timescale for the magnetic energy to adjust 
toward its expectation value, i^iss is the timescale for 
magnetic dissipation, and tcooi is the radiation cooling 
time. The magnetic fleld energy matches the equilib- 
rium value when ^growth = idiss, so these two timescales 
may be equated. Non-dimensionalizing the energies in 
terms of their equilibrium values and the time in terms 
of idiss, these equations become 

'^^^ :i?(i)£]^-£B (23) 



dr 
d£R 
dr 



Ebo idiss c 
-^B— -; OR. 



(24) 



Ero icool 

The timescale ratio idiss/^cooi may also depend on pres- 
sure; we designate its logarithmic derivative with respect 



to pressure by — s (the minus s ign makes its deflnition 
consis tent with the notation in iKrolik. Hirose. fc BlaesI 
(I2OOI ). 

That assumption puts the model equations in the 
form: 



d£B 

d£R 
dr 



R{t)£^ — Eb 
Ebo 



E 



^diss,0 pi — s 



RO 



tr. 



□1,0 



(25) 
(26) 



Once again we can use the requirement that equilib- 
rium values be consistent with zero time-dependence to 
constrain these parameters; this time we conclude that 
idiss,o/^cooi,o = Ebo/Ero, leaving 



- Rit)E'^ - Eb 
dER, Ebo [c c 



(27) 
(28) 



In other words, the timescale for the radiation energy to 
equilibrate is ~ Erq/Ebo times the timescale on which 
the magnetic energy can fluctuate; the magnetic energy 
must be dissipated many times over in order to produce 
a response in the radiation energy because Erq/Ebo is 
so large. For context, we recall that in the simulations 
described here, the mean Er/Eb — 30 while the mean 
thermal timescale was ~ 15 orbits. If this toy-model 
describes the simulations, i.e., if heating in the simula- 
tions is primarily due to magnetic dissipation and cooling 
is due to radiation losses, these numbers would suggest 
that in the simulations tdiss should be defined as the ratio 
of box-integrated magnetic energy to magnetic dissipa- 
tion rate, and its typical value should be ~ 0.5 orbits. In 
fact, consistent with this description of energy balance, 
the mean value of idiss defined in this way was 0.53 orbits 
in simulation 1112a, and 0.55 orbits in simulation 1126b. 

In Fig. [2TI we show a sample solution of these equa- 
tions, in this case with n = s = 0, Ebo/Erq = 0.033, 
and initial conditions Er = Eb = ^- Note that setting 
n = means that there is no explicit dependence of mag- 
netic amplification on the radiation pressure. The ran- 
dom function R{t) was constructed by requiring its power 
spectrum to match approximately the power spectrum 
of magnetic energy fluctuations seen in simulation 1112a 
(i.e., oc f~^'^^ at frequencies / < 2.6/tcooi and oc /-3.65 
for / > 2.6/icooi), but choosing the phases of its Fourier 
amplitudes randomly with uniform probability distribu- 
tion in the interval [0,27r). To make the comparision as 
direct as possible, the duration of the integration was set 
to 40icooi- Just as in the simulations, the radiation en- 
ergy fluctuates by factors of several over this timespan, 
but shows no clear trend. 

Guided by a traditional view of the a model, one might 
object that the absence of thermal instability here is due 
to the lack of any explicit dependence of the magnetic 
energy on the radiation energy. However, closer study 
of the differential equations suggests that if icooi is short 
compared to the timescale of the large-amplitude turbu- 
lence fluctuations, approximate radiative equilibrium is 
enforced on the shorter timescale. It would then follow 
that Eb cx: fjj"" independent of n. This is exactly what 
we see, for these parameters and all other combinations 
of n and s we have tried, subject only to the constraint 
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Fig. 21. — Radiation energy (black curve) and magnetic energy 
(gray curve) as functions of time in units of tcooli computed on 
the basis of the toy-model described in the text. 
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Fig. 22. — Correlation plot between Sn and £b- 

that n < 1 — s (as discussed in iKrolik. Hirose. fc Blaed 
()2007! ). genuine instability can be found when this limit is 
violated). The flat power spectrum at frequencies below 
l/^cooi guarantees that there is substantial fluctuation 
power on timescales longer than the radiative equilibra- 
tion timescale, and the result is as predicted, as shown 
by Fig. [22l A linear least-squares flt in the logarithms 
of these two quantities then yields a best-flt relation in 
which d\og£B/d^og£ji ~ 1 — s; in the example shown, in 
which s = 0, the best-flt logarithmic derivative was 1.01. 
Thus, magnetic energy and pressure do scale together as 
suggested by the dimensional analysis underlying the a 
model, but it is not necessarily because the pressure di- 
rectly forces the magnetic energy; it is instead the other 
way around: the pressure is the result of magnetic dissi- 
pative heating regulated by photon losses. 

We now see the data shown in Fig. [S] in a different 
light. Although it is true that, compared at the same 
time, Eb oc E^''^ in that simulation, this correlation 
is created by the dependence of the radiation energy 
on the magnetic energy, and not by the radiation pres- 
sure regulating the magnetic fleld strength. Moreover, 
if our simple model correctly describes the thermal bal- 
ance, the slope of this correlation should be 1 — s, so 
that tcooi/idiss (X E'^'^^. We have already seen that 
tcooi c!c Ej^^^; a logarithmic least-squares linear flt to 
the data of 1112a gives ^diss oc supporting the toy- 

model. 

An important corollary of these results is that the in- 



trinsic dependence of magnetic energy on radiation pres- 
sure parameterized by n cannot be inferred from the slope 
of the correlation between magnetic energy and radia- 
tion energy; this slope is entirely independent of it. The 
most that one can do is to place an upper bound on how 
strongly pressure influences magnetic saturation by using 
the observed thermal stability to argue that n < 1 — s. 

To close this section, we make the qualitative remark 
that the sensitivity of the cooling time to pressure may 
influence the magni tude of fluctuations in the energy con- 
tent. As argued in lKrolik. Hirose, fc BlaesI (pb07,l . pro- 
vided idiss is not too strong a function of box energy 
content, s must be relatively large in magnitude and 
negative when gas pressure dominates radiation pressure 
because an increasing ratio of radiation to gas pressure 
makes diffusive cooling much more rapid. If so, the cool- 
ing time rapidly shortens when the radiation energy in- 
creases, placing a tight cap on fluctuations in the energy 
content of the box. When the radiation to gas pres- 
sure ratio is order unity or greater, the scaling of <cooi 
with total energy should be rather slower (a prediction 
consistent with our results if idiss is likewise not a strong 
function of Ef;). In this case, the cooling time is compar- 
atively insensitive to energy content and exercises looser 
control on energy fluctuations, permitting the turbulence 
fluctuations to drive larger amplitude fluctuations in the 
total pressure. 

5. CONCLUSIONS AND SUMMARY 

In this paper we have presented the results of two sim- 
ulations that each followed the evolution of a vertically- 
stratifled shearing box with the same surface density and 
orbital frequency for 600 orbits, or ~ 40 thermal times. 
The initial conditions of the two simulations differed only 
in being given different realizations of the small ampli- 
tude noise that was imposed on their otherwise smooth 
initial state. Because the magneto-rotational instability 
generates MHD turbulence, these systems are chaotic; 
differing small amplitude noise therefore leads to order 
unity contrasts in their subsequent evolution. The two 
simulations should thus be viewed as two different real- 
izations of the many evolutions possible for these param- 
eters. Their qualitative aspects are, nonetheless, similar. 

For example, just as we have found in previ- 
ous simulations studying shearing boxes with p,. ~ 
0.2pp (Hirose, Krolik, & Stone 2006) and Pr ~ Pg 
(jKrolik, Hirose, fc BlaesI I2007f ). most of the disk mass is 
found near the midplane, where the magnetic fleld en- 
ergy is ^ 10^^ ti mes the com bined gas and radiation 
pressure (see also iTurneil |2004| ). Likewise, in all cases 
the upper layers of the disk are magnetically-dominated, 
and the photosphere lies within this region. Unlike the 
earlier simulations, in these two the box-averaged radia- 
tion pressure is ~ 10 times greater than the gas pressure 
at all times. 

Most importantly, in both cases, the energy con- 
tent of the box undergoes order unity fluctuations over 
timescales of many tens of orbits, but these fluctuations 
have no long-term tren d. This result directly chall enges 
the prediction made in IShakura fc SunvaevI ()19760 that 
radiation pressure-dominated disk segments should be 
the rmally unstable. Not e , how ever, that the prediction 
by iLightman fc EardlevI ()1974l ) of inflow instability in 
radiation-dominated disks remains to be investigated. 
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In our simulations, the time- averaged dissipation rate 
in the disk body is roughly equal to the characteris- 
tic value cfi^/Kcs, the rate at whic h diffusive radiation 
flux m aintains hydrostatic balance (jShakura fc SunvaevI 
[1976"). However, in some places (|z| ~ H), the local 
dissipation rate exceeds the characteristic value by more 
than 30%. If this heat went into radiation flux, hydro- 
static balance would be disrupted. We find that the ex- 
cess is exactly compensated by radiative advection as- 
sociated with an acoustic breathing mode. At the same 
time, mechanical work done on the fluid by the shearing 
boundaries in excess of the local dissipation rate is trans- 
ported outward by Poynting flux and radiation pressure 
work associated with the breathing mode. Both the dis- 
sipated energy carried by radiation advection and (non- 
dissipated) energy carried by Poynting flux and radia- 
tion pressure work are finally deposited and dissipated 
at \z\ ~ 2H, and from there all the way through the rest 
of the structure radiation diffusion overwhelms all other 
energy fluxes. 

To explain the thermal stability of radiation- 
dominated disks, we argue that the comparability of 
stress and pressure inferred from dimensional analysis 
(which underlies the a-model, and the prediction of in- 
stability) is due to dissipation of magnetic turbulence 
(which produces the stress) providing the heat that is 
then transformed into radiation pressure. Consequently, 
fluctuations in magnetic energy drive fluctuations in 
pressure, and not (as has been commonly assumed) the 
other way around. Our claim is supported by two lines 
of evidence: First, crosscorrelation analysis of simula- 
tion data demonstrates that magnetic fluctuations lead 
radiation energy fluctuations by 5-15 orbits, a little less 
than a thermal time. Pressure fluctuations cannot, then, 
drive magnetic fluctuations. Second, we constructed a 
simple model realizing this picture, and this model re- 
produces two other important features observed in the 
simulation data: The system undergoes thermal fluctua- 
tions closely resembling in amplitude and timescale those 
seen in the simulations. In addition, although no corre- 
lation between magnetic energy and radiation pressure is 
built into the model, thermal balance automatically cre- 
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ates one after the fact. Thus, correlations between stress 
and pressure are due to the dissipation of magnetic en- 
ergy supplying thermal energy, not to the pressure defln- 
ing a characteristic scale for the stress. 

A logical consequence of this point of view, in which 
stress determines pressure, rather than the other way 
around, is that the fundamental independent variables 
are surface density and orbital frequency. That the or- 
bital frequency is independent of disk parameters is ob- 
vious, so long as the disk mass is small compared to the 
central mass. So long as the inflow timescale is long 
compared to the thermal timescale, the surface density 
must likewise be regarded as an independent parame- 
ter with respect to thermal and dynamical fluctuations. 
These two independent parameters, through the inter- 
twined and nonlinear processes of MHD instability, tap- 
ping the energy reservoir of orbital shear, magnetic dis- 
sipation, thermal radiation, and radiative diffusion, with 
all of these taking place under conditions of vertical (as 
well as radial) gravitational confinement, combine to de- 
termine the magnetic field strength, both its mean satu- 
ration level and its fluctuations. Orbital shear flxes the 
stress exerted by this field, while the turbulent cascade 
sets the dissipation rate. These two are not entirely sep- 
arate, of course, as the time-averaged accretion rate is 
equivalent to either one. The pressure follows from the 
heating rate, as regulated by photon diffusion, and, in 
turn, closes the loop by determining the disk thickness. 
Despite all these complications, at bottom, everything 
is still determined by only two variables, surface density 
and orbital frequency. 
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APPENDIX 

NUMERICAL METHODS FOR SOLVING THE MATTER-RADIATION ENERGY EXCHANGE 

EQUATIONS 

To solve the equations for matter energy density e (eqn. [S]) and radiation energy density E (eqn. [S]), we follow our 
usual system of operator-splitting. Because energy removed from one of these reservoirs can go directly into the other, 
we treat the operator-split segments in pairs. 

Free-free absorption 

First we group together the work done by matter on radiation and the energy exchanged by free-free absorption and 
emission: 

de 

— = -{y.v)p~{A^B-cE)Rlp (Al) 

— = -Vv : P + {inB - cE)R%p. (A2) 

To solve these equations implicitly, we follow the method described in § 4.3 in [Turner fc Stoii3 (|2001[ ). where the 
following quartic for x = e"+^ is derived (the superscript is the time-step index): 

x'^ + C2X + ci=Q (A3) 



where 



_ (l + a3-fa2)e"+a2ig" , . 

ci = 7T— — ^ , (A4) 

_ (1 + 04X1 + 03 + 02) .... 

C2 = ZT- T , ( A5) 



01(1 -I- 03) 

m(7-1) 



4 



ai=44aB [^ ^'-i^pn" j (A6) 
a2 = cK^At, (A7) 

«3^ ^ J, J, ^t, (A8) 

a4 = {-/-l){V ■v)At. (A9) 

Here TZ and ctb are the gas constant and the Stefan-Boltzmann constant, respectively. Once the quartic (jA3p is solved 
for X = e"+i, E"+^ is computed as {E" + e" - {I + a4)x)/{l + 03). 

Hereafter we assume that (l-t-aa) = {E"+'^ + (yv : P)"+'^ At) / E"-+'^ > and (l+Ui) = (e"+i-h(V-t>)p"+iAi)/e"+i > 
0. This assumption is valid whenever At is small enough for the changes in the energy densities by photon damping 
and gas compression during the time-step to be smaller than the energies themselves. Then ci < and C2 > because 
oi > and 02 > 0. The signs of ci and C2 are crucial in the following.^ 



Solving the quartic lfA3\} by Ferrari's formula 
The quartic (|A3[) can be solved by Ferrari's formula. First we rewrite the equation as 

X* = -C2X-Ci. (AlO) 

Adding {tx^ + t^/4) to both sides leads to a new quartic: 

tV . e 



a;2 + -) = -cix-cx^tx^ (All) 



If we could find a value t such that 



we could rewrite equation (|All[) as 



C2X - ci + te"^ + J = [\rtx - , (A12) 



2/ V " 2^/^ 
An equation in this form can be solved readily, giving the four roots 



±\/t+ V-<T(2c2)/\/i ±Vi- V-iT(2c2)/v^ 
2^1,2 = , X3,4 = ^—^ . (A14) 

^ In our actual simulations, we exclude the gas compression term and solve for its effect separately, which means 04 = 0, but (l-l-as) > 
is always true. If the photon damping term is also treated separately (03 = 0), it is always guaranteed that ci < and C2 > 0. 
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Since ci < and C2 > as discussed above, there is only one positive solution for x: 



-Vi+J-t+{2c2)/^t 

^- . (A15) 



2 

If we can assume that t > 0, this solution is real. We next prove that this is a good assumption. 

Solving the cubic SAISH by Cardano 's formula 
Equation IA12I can be rewritten as a cubic for t: 

+ 3pt + q^ 0, (A16) 

where p = — 4ci/3 and q = —(? . Here p > and g < because ci < and ci > 0. The single positive solution of this 
cubic can be obtained via Cardano's formula as follows: 

t ^ ^/a^ + .^/a^, (A17) 

where 

a± = . (A18) 

Note that i > 0, as required. 

Compton scattering 

Next consider the terms describing Compton scattering. We rewrite them in the form 

dE ( 4crps \ ^ . , ( E'^ 



.E\(n- l)Aimpe - fcBP - (A19) 
ot \ ficm-cmpC J \ \a J I 

Here Cos, ?7ip, and nic are the Thomson cross section, the mass of the proton, and the mass of electron, respectively. We 
assume that the mean molecular weights fi and fi^ are 0.61 and 1.2 respectively. We solve these equations implicitly: 

' ^i?"+M(7-lVmpe"+i~A:Bp"(^) (A21) 



e 



= 0. (A22) 



At At 

Eliminating e"+^ from this pair of equations, we find a nonlinear equation for x = (i?"+^/_E")^/'': 



x''l^x^ + Ax+i^B-l-—jj-B = 0, (A23) 

where A = ( fcBp"(-E"/a)^/'^)/((7 - l)fimpE'^) and B = (//eTOec)/(4^(7 - l)aesE"At). We solve the nonlinear 
equation IA23I by the Newton-Raphson method. That gives us the energies for the next step, £'"+^ = and 
e"+i = (£"» + e") - 



